library(rgrass7)
## Loading required package: XML
## GRASS GIS interface loaded with GRASS version: GRASS 7.8.2 (2019)
## and location: rdom
gisdbase <- 'grass-data-test' #Base de datos de GRASS GIS
wd <- getwd() #Directorio de trabajo
wd
## [1] "/home/edelte/unidad-0-asignacion-99-mi-manuscrito-edeltejeda"
loc <- initGRASS(gisBase = "/usr/lib/grass78/",
home = wd,
gisDbase = paste(wd, gisdbase, sep = '/'),
location = 'c_ozama',
mapset = "PERMANENT",
override = TRUE)
gmeta()
## gisdbase /home/edelte/unidad-0-asignacion-99-mi-manuscrito-edeltejeda/grass-data-test
## location c_ozama
## mapset PERMANENT
## rows 1
## columns 1
## north 1
## south 0
## west 0
## east 1
## nsres 1
## ewres 1
## projection NA
execGRASS(
'g.list',
flags = 't',
parameters = list(
type = c('raster', 'vector')
)
)library(raster)
## Loading required package: sp
rutadem <- 'data/dem.tif'
rawextent <- extent(raster(rutadem))
rawextent
## class : Extent
## xmin : 314490.3
## xmax : 460848.8
## ymin : 2008160
## ymax : 2103177
devtools::source_url('https://raw.githubusercontent.com/geofis/rgrass/master/integerextent.R')
## SHA-1 hash of file is 65b57fe7cb20cd1976572cd0a7e98def9e95d83c
devtools::source_url('https://raw.githubusercontent.com/geofis/rgrass/master/xyvector.R')
## SHA-1 hash of file is 89fa5ae436d9a7d7a0c799b789c560eb5e421cfd
newextent <- intext(e = rawextent, r = 90, type = 'inner')
newextent
## class : Extent
## xmin : 314550
## xmax : 460800
## ymin : 2008170
## ymax : 2103120
gdalUtils::gdalwarp(
srcfile = 'data/dem.tif',
dstfile = 'data/demint.tif',
te = xyvector(newextent),
tr = c(90,90),
r = 'bilinear',
overwrite = T
)
## Registered S3 method overwritten by 'R.oo':
## method from
## throw.default R.methodsS3
## NULLrutademint <- 'data/demint.tif'
execGRASS(
"g.proj",
flags = c('t','c'),
georef=rutademint)
gmeta()
execGRASS(
"r.in.gdal",
flags='overwrite',
parameters=list(
input=rutademint,
output="demint"
)
)
execGRASS(
"g.region",
parameters=list(
raster = "demint",
align = "demint"
)
)gmeta()
## gisdbase /home/edelte/unidad-0-asignacion-99-mi-manuscrito-edeltejeda/grass-data-test
## location c_ozama
## mapset PERMANENT
## rows 1055
## columns 1625
## north 2103120
## south 2008170
## west 314550
## east 460800
## nsres 90
## ewres 90
## projection +proj=utm +no_defs +zone=19 +a=6378137 +rf=298.257223563
## +towgs84=0.000,0.000,0.000 +type=crs +to_meter=1
execGRASS(
'g.list',
flags = 't',
parameters = list(
type = c('raster', 'vector')
)
)
## raster/demintexecGRASS(
"r.stream.extract",
flags = c('overwrite','quiet'),
parameters = list(
elevation = 'demint',
threshold = 80,
stream_raster = 'stream-de-rstr',
stream_vector = 'stream_de_rstr'
)
)
execGRASS(
'g.list',
flags = 't',
parameters = list(
type = c('raster', 'vector')
)
)
## raster/demint
## raster/stream-de-rstr
## vector/stream_de_rstrlibrary(sp)
use_sp()
library(mapview)
netw <- spTransform(
readVECT('stream_de_rstr'),
CRSobj = CRS("+init=epsg:4326"))mapview(netw, col.regions = 'blue', legend = FALSE)
## Warning in CPL_crs_from_input(x): GDAL Message 1: +init=epsg:XXXX syntax is
## deprecated. It might return a CRS with a non-EPSG compliant axis order.
## Warning: sf layer has inconsistent datum (+proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs).
## Need '+proj=longlat +datum=WGS84'source('my-trans.R')
outlet <- as.integer(my_trans(c(-69.88117,18.47503)))r.basinpref <- 'rbasin_ozama'
execGRASS(
"r.basin",
flags = 'overwrite',
parameters = list(
map = 'demint',
prefix = pref,
coordinates = outlet,
threshold = 80,
dir = 'salidas-rbasin/ozama'
)
)execGRASS(
'g.list',
flags = 't',
parameters = list(
type = c('raster', 'vector')
)
)
## raster/demint
## raster/rbasin_ozama_demint_accumulation
## raster/rbasin_ozama_demint_aspect
## raster/rbasin_ozama_demint_dist2out
## raster/rbasin_ozama_demint_drainage
## raster/rbasin_ozama_demint_hack
## raster/rbasin_ozama_demint_hillslope_distance
## raster/rbasin_ozama_demint_horton
## raster/rbasin_ozama_demint_shreve
## raster/rbasin_ozama_demint_slope
## raster/rbasin_ozama_demint_strahler
## raster/stream-de-rstr
## vector/rbasin_ozama_demint_basin
## vector/rbasin_ozama_demint_mainchannel
## vector/rbasin_ozama_demint_network
## vector/rbasin_ozama_demint_outlet
## vector/rbasin_ozama_demint_outlet_snap
## vector/stream_de_rstrSi
r.basinarrojara error (sólo en el caso de error, no en caso de advertencia), ejecutar este bloque para borrar las salidas anteriores y reejecutar elr.basin:
execGRASS(
"g.remove",
flags = 'f',
parameters = list(
type = c('raster','vector'),
pattern = paste0(pref, '*')
)
)rbnetw <- spTransform(
readVECT('rbasin_ozama_demint_network'),
CRSobj = CRS("+init=epsg:4326"))
rbnetw
rbmain <- spTransform(
readVECT('rbasin_ozama_demint_mainchannel'),
CRSobj = CRS("+init=epsg:4326"))
rbmain
rbbasin <- spTransform(
readVECT('rbasin_ozama_demint_basin'),
CRSobj = CRS("+init=epsg:4326"))
rbbasinlibrary(leaflet)
leaflet() %>%
addProviderTiles(providers$Stamen.Terrain, group = 'terrain') %>%
addPolylines(data = rbnetw, weight = 3, opacity = 0.7) %>%
addPolylines(data = rbmain, weight = 3, opacity = 0.7, color = 'red') %>%
addPolygons(data = rbbasin) %>%
leafem::addHomeButton(extent(rbbasin), 'Ver cuenca')library(readr)
rbozamapar1 <- read_csv("salidas-rbasin/ozama/rbasin_ozama_demint_parametersT.csv")
## Parsed with column specification:
## cols(
## .default = col_double(),
## Rectangle_containing_basin_N_W = col_character(),
## Rectangle_containing_basin_S_E = col_character()
## )
## See spec(...) for full column specifications.
rbozamapar1 %>% tibble::as_tibble()
## # A tibble: 1 x 34
## x y Easting_Centroi… Northing_Centro… Rectangle_conta…
## <dbl> <dbl> <dbl> <dbl> <chr>
## 1 406935 2.04e6 413955 2069865 ('367560', '209…
## # … with 29 more variables: Rectangle_containing_basin_S_E <chr>,
## # Area_of_basin_km2 <dbl>, Perimeter_of_basin_km <dbl>,
## # Max_Elevation <dbl>, Min_Elevation <dbl>, Elevation_Difference <dbl>,
## # Mean_Elevation <dbl>, Mean_Slope <dbl>,
## # Length_of_Directing_Vector_km <dbl>,
## # Prevalent_Orientation_deg_from_north_ccw <dbl>,
## # Compactness_Coefficient <dbl>, Circularity_Ratio <dbl>,
## # Topological_Diameter <dbl>, Elongation_Ratio <dbl>,
## # Shape_Factor <dbl>, Concentration_Time_hr <dbl>,
## # Length_of_Mainchannel_km <dbl>,
## # Mean_slope_of_mainchannel_percent <dbl>,
## # Mean_hillslope_length_m <dbl>, Magnitudo <dbl>,
## # Max_order_Strahler <dbl>, Number_of_streams <dbl>,
## # Total_Stream_Length_km <dbl>, First_order_stream_frequency <dbl>,
## # Drainage_Density_km_over_km2 <dbl>, Bifurcation_Ratio_Horton <dbl>,
## # Length_Ratio_Horton <dbl>, Area_ratio_Horton <dbl>,
## # Slope_ratio_Horton <dbl>
rbozamapar2 <- read_csv(
"salidas-rbasin/ozama/rbasin_ozama_demint_parameters.csv",
skip=2, col_names = c('Parameter', 'Value'))
## Parsed with column specification:
## cols(
## Parameter = col_character(),
## Value = col_character()
## )
rbozamapar2 %>% print(n=Inf)
## # A tibble: 32 x 2
## Parameter Value
## <chr> <chr>
## 1 Easting Centroid of basin 413955.00
## 2 Northing Centroid of basin 2069865.00
## 3 Rectangle containing basin N-W ('367560', '20983…
## 4 Rectangle containing basin S-E ('455220', '20399…
## 5 Area of basin [km^2] 3408.388875
## 6 Perimeter of basin [km] 358.481006282309
## 7 Max Elevation [m s.l.m.] 913.3936
## 8 Min Elevation [m s.l.m.] -3.278881
## 9 Elevation Difference [m] 916.672481
## 10 Mean Elevation 106.245
## 11 Mean Slope 3.47
## 12 Length of Directing Vector [km] 27.788971913332816
## 13 Prevalent Orientation [degree from north, counterclo… 1.3166010130150796
## 14 Compactness Coefficient 5.441724127967625
## 15 Circularity Ratio 0.333293391887047…
## 16 Topological Diameter 168.0
## 17 Elongation Ratio 0.5008280978139896
## 18 Shape Factor 25.91243324088714
## 19 Concentration Time (Giandotti, 1934) [hr] 17.787167537498174
## 20 Length of Mainchannel [km] 131.534883016
## 21 Mean slope of mainchannel [percent] 1.1417035909306994
## 22 Mean hillslope length [m] 4962.319
## 23 Magnitudo 1096.0
## 24 Max order (Strahler) 7
## 25 Number of streams 1615
## 26 Total Stream Length [km] 3117.2892
## 27 First order stream frequency 0.321559552091895…
## 28 Drainage Density [km/km^2] 0.9145931741723573
## 29 Bifurcation Ratio (Horton) 3.3679
## 30 Length Ratio (Horton) 1.6319
## 31 Area ratio (Horton) 3.6528
## 32 Slope ratio (Horton) 1.3764unlink_.gislock()